//      dna_processing.cc
//      
//      Copyright 2012 tqlong <tqlong@espada>
//      
//      This program is free software; you can redistribute it and/or modify
//      it under the terms of the GNU General Public License as published by
//      the Free Software Foundation; either version 2 of the License, or
//      (at your option) any later version.
//      
//      This program is distributed in the hope that it will be useful,
//      but WITHOUT ANY WARRANTY; without even the implied warranty of
//      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//      GNU General Public License for more details.
//      
//      You should have received a copy of the GNU General Public License
//      along with this program; if not, write to the Free Software
//      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
//      MA 02110-1301, USA.

#include <cctype>
#include "stl_helper.h"
#include "dna_processing.h"

int Sequences::max_length() const
{
  return 0;
}

void Sequences::fasta_read(const std::string& fastaFile)
{
  FileAccessRead(file, fastaFile);
  const int BUF_SIZE = 1024*1024;
  char buf[BUF_SIZE];
  int n, s = 0; // buf[s] <-- current character
  std::string a(1024*1024, ' '), seq(1024*1024, ' ');
  int state = 0;

  //MSG("Read 1M");
  file.read(buf, BUF_SIZE);
  n = file.gcount();
  int pos = 0, pos_anno = 0;
  this->n = this->n_gc = 0;
  while (n > s) {
    //MSG("state=" << state << " " << n << " " << s << " seq=" << seq << " a=" << a);
    switch (state) {
    case 0: if (buf[s] == '>') state = 1; break;
    case 1: 
      if (buf[s] != '\n' && buf[s] != '\r') {
	a[pos_anno++] = buf[s];
      }
      else {
	state = 2;
	pos = 0;
      }
      break;
    case 2:
      if (buf[s] == 'A' || buf[s] == 'C' || buf[s] == 'G' || buf[s] == 'T' ||
	  buf[s] == 'a' || buf[s] == 'c' || buf[s] == 'g' || buf[s] == 't'
	  ) {
	this->n++;
	if ( buf[s] == 'C' || buf[s] == 'G' || buf[s] == 'c' || buf[s] == 'g' ) this->n_gc++;
	seq[pos++] = (char)toupper(buf[s]);
	if (pos >= sz(seq)) {
	  //MSG("resize " << pos << " " << sz(seq));
	  seq.resize(2*sz(seq), ' ');	  
	}
      }
      else if (buf[s] == '>') { // next sequence detected
	this->pb(seq.substr(0, pos));
	this->annotation_.pb(a.substr(0, pos_anno));
	pos_anno = 0;
	state = 1;
      }
      break;
    default:
      FATAL_ERROR("State not recognized");
    }

    s++;
    if (s >= n) {
      //MSG("Read 1M");
      file.read(buf, BUF_SIZE);
      n = file.gcount();
      s = 0;
    }
  }

  if (pos > 0) { // push_back the last one
    this->pb(seq.substr(0,pos));
    this->annotation_.pb(a.substr(0,pos_anno));
  }
}

std::ostream& Sequences::print (std::ostream&o, bool isFile) const
{
  for (uint i = 0; i < this->size(); i++) {
    o << ">" << annotation_[i] << ENDL;
    if (!isFile) 
      o << this->at(i) << ENDL;
    else {
      for (uint j = 0; j < this->at(i).size(); j += 80) {
	int s = 80;
	if (j + s > this->at(i).size()) s = this->at(i).size() - j;
	o << this->at(i).substr(j, s) << ENDL;
      }
    }
  }
  return o;
}

void Sequences::split(std::ostream&o, int length)
{
  for (int i = 0; i < sz(*this); i++) {
    int count = 0;
    const std::string& seq = this->at(i);

    for (int j = 0; j <= sz(seq)-length; j += length) {
      o << ">" << i << " " << count++ << ENDL;
      o << seq.substr(j, length) << ENDL;
    }
  }
}
